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Abstract 

In this paper un-binned statistical tools for analyzing the cosmic ray energy spec- 
trum are developed and illustrated with a simulated data set. The methods are 
designed to extract accurate and precise model parameter estimators in the pres- 
ence of statistical and systematic energy errors. Two robust methods are used to 
test for the presence of flux suppression at the highest energies: the Tail- Power 
statistic and a likelihood ratio test. Both tests give evidence of flux suppression in 
the simulated data. The tools presented can be generalized for use on any astro- 
physical data set where the power-law assumption is relevant and can be used to 
aid observational design. 

Key words: cosmic ray spectrum, power-law, CRPropa, TP-statistic, flux 
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1 Introduction 



The observation of suppression in the flux of the highest energy cosmic rays 
(CRs) has been of central interest to astro-particle physics since the predic- 
tion of the GZK-effect |6|ll7j in 1966. Most recently both the Auger [15] and 
the HiRes[T] detectors have released results favoring the observation of flux 
suppression at a 6cr and 5a level of confidence, respectively. 

With this in mind, we describe a set of statistical tools designed to extract 
the most accurate and precise information concerning the flux of the highest 
energy cosmic rays. By binning the data we can only lose information [5J (see 
CA|) and therefore our statistical tools use an un-binned maximum likelihood 

^ Corresponding author, E-mail: jhagueOunm.edu 
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approach [IGMTTP] to answer two related statistical questions: Is there flux 
suppression at the highest energies? and, if yes, What are the characteristic 
cut-off energy and shape parameters? 

In detail we first generate a toy data set using the CRPropa package[2], as in 
§2.21 We then fit this simulated data to the three models described in §2.31 
The un-binned maximum likelihood fit is outlined in §3.11 and methods for 
incorporating systematic and statistical energy errors are described in §3.21 
and §3.31 respectively. In ^we describe several statistical tools for hypothesis 
testing: the Kolmogorov-Smirnov test, the tail power statistic p!2l[71[T5] . and a 
likelihood ratio test[S]. 

Though we cast our discussion in terms of cosmic ray energies, it is worth 
noting that these tools can be applied to any astrophysical data set where de- 
viations from the power-law hypothesis are relevant, e.g. the galaxy correlation 
function |18j or gamma ray astronomy |13j. 



2 CRPropa Data Set and Models 

2.1 Input from the HiRes and Auger Observatories 

Both the HiRes [1] and Auger [15] observatories have reported spectra and fit 
parameters for various power-law models. The collaborations use binned fitting 
methods. They fit the spectrum over many orders of magnitude in energy 
but we summarize here the model parameter^ relevant only to the highest 
energies. The best fit double power-law parameters reported by HiRes [IJ are 
7 = 2.81±0.03(stat)±0.02(sys), = 10^-^5±°-°^(stat) and 5 = 5.1 ±0.7(stat). 
For the same model Auger [T5] reports 7 = 2.62 ± 0.03(stat)±0.02(sys), E]^ = 
lO^'^(fixed) and 5 = 4.14±0.42(stat). Fitting to the Fermi power-law Auger [T5] 
finds 7 = 2.56±0.06(stat), Ei = 10i-^^^°-°^(stat) and Wc = 0.16 ± 0.04(stat). 

2.2 A Toy CR Data Set 

To illustrate the methods in this note we use un-binned proton primary cos- 
mic ray, CR, arrival energies (in EeV= lO^^eV) as simulated by the pack- 
age CRPropa [2j with input spectral index 71N = 2.6, -Emin = 10 EeV and 
-E'max = 2000 EeV. We draw 5 x 10^ events to act as a toy data set from a 
modern CR detector. 



^ See ^2.31 and Tabled] for the definition of these parameters. 
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Fig. 1. The differential flux as simulated by 5 x 10^ events from the CRPropa toy 
set with parameters 71N = 2.6 and -Emax = 2000 EeV (see §2.2p . The p.d.f. of the 
best fit double power-laws reported by HiRes [Ij and Auger [15j are the dashed lines. 



The CRPropa toy data set is similar size and shape to the flux reported by 
these observatories but the results of this study do not, otherwise, reflect any 
information about any physical data set. The probability distribution function 
(p.d.f.) of the best fit double power-laws reported by HiRes[lJ and Auger [T3] 
are shown in Fig{T] along with the CRPropa toy data. 

The CRPropa propagation simulation is implemented by first generating pro- 
ton CR primaries with initial energies according to a power-law "at the source," 
propagating them through a simulated Universe and then observing the final 
energy. The spacial extent of the sources is simulated as a uniform distribution 
of discrete sources on a grid with 10 Mpc steps extending to a distance of 4.07 
Gpc, (from redshift z = 0.0 to z = 2.73). Nuclei traveling over many mega- 
parsecs from these sources will suffer significant energy loss in an expanding 
Universe filled with the cosmic microwave background, CMB, radiation. As 
a result, the highest energy flux is much less than one would expect from a 
power-law alone. This suppression is known as the GZK-effect [61IT7] . 



2.3 Power-Law Models 



The fundamental probability distribution function governing the pure power- 
law assumption, denoted /p, is shown in Table [U /p = (7 — l)E^in^~~^E~'^ . 
The parameter 7 is referred to as the spectral index. Here the sub-scripted-P 
stands for Pure-power-law. 

For the highest energy CRs, the interesting observation would be to confirm 
or deny deviation from the power-law form at the highest magnitudes, i.e. the 
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Fig. 2. The differential flux as simulated by 5 x 10^ events from the CRPropa toy set 
with parameters 7in = 2.6 and -Emax = 2000 EeV. The best fit models are described 
in g231 



GZK-cutoff. We therefore study two toy models that mimic a pure power-law 
for lower energies but exhibit flux suppression above a given energy. The first 
is a double power-law (DP) with two spectral indexes, 7 below Ei, ("b" for 
bend or break) and 5 > 7 above. The point at which this p.d.f. reaches half 
the value it would have if the pure power-law continued above i?b is given by 
see [3] for a discussion of this quantity. Both HiResp^J and 
Auger [15] have analyzed their data using this model. 

We also study a toy p.d.f. where the cut-off is a "Fermi-like" Power-law 
(FP) [T5|7] . The advantage of fitting with this toy model is that the location 
parameter Ei is a parameter in the fit. 

All three p.d.f. 's are normalized on the interval [E^^^, 00), i.e. 
()m = ISL,jM{t)dt = 1 for each of the models M e {P,DP,FP}. The first 
element of the parameter vector 9i = -Emin is fixed for the fit (see ^ and 
then varied to estimate the stability (see §4.ip . Thus the power-law has one 
free parameter and the other models have three; low energy spectral index, 
location of cut-off and "steepness" of cut-off. 



3 Fitting the Data 



We take an un-binned maximum log-likelihood approach to estimating the 
best-fit parameters of each model. The method constructed here is designed to 
extract the maximum possible statistical information about these parameters. 
For the ideal detector we assume that the observed energies are known with 
infinite precision. 
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The model designation (Model = Pure power-law, Double Power-law or Fermi 
Power-law), number of free parameters, normalization, and form of the function 
used to fit the simulated fluxes used in this study. 



3. 1 Ideal Detector 

We find estimates of the parameters in each model by maximizing, 

N 

CM{e) = Y.\n[fM{Ej)], (1) 

i=l 

where the sum is carried out over the event energies and 9i = -Emin is fixed. 
The global maximum of this function Cu{(^) determines the best parameter 
estimates, 6. The the function is maximized using Minuit[in] with the MIGrad 
option. 

To determine the one degree of freedom error estimate [I6] for a parameter 
we vary the parameter (with the others fixed at 6) until — 2A£m = 1- The 
two degrees of freedom error estimates [T6] are determined by varying two 
parameters with the other fixed and choosing the contour such that —2ACm > 
2.3. For the toy data set, we plot these contours and the asymmetric one degree 
of freedom error estimates in ^|Cl Fig JC.3l and IC.4[ 

3.2 Systematic Energy Error 

The errors on the observed energy Eobs of an event from a real CR detector 
are considerable and must be included in any realistic analysis of a spectrum. 
For our purposes, these errors take the two canonical forms; statistical and 
systematic, i.e. E^bs ± c^stat ± (^sys- 

The systematic errors energy errors of a CR detector reflect the uncertainties 
in the absolute calibration of the detector. At the highest energies the system- 
atics are the dominant contribution to the overall uncertainty of an event's 
energy. For example, the two fluorescence detectors Auger 115] and Hires[lJ re- 
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port uncertainties of 22% and 17% respect ivelj{J] The shift in energy due to the 
systematic error can be asymmetric, i.e. cr^g 7^ o"^^, and energy dependent, 
see Eq([2]), but it effects every event at a given energy the same way; a shift 
up or down. For the Monte-Carlo (MC) data sets we model the systematic 
detector energy errors using: 

^^=Pi+P2HE). (2) 

Here we choose symmetric systematically-shifted energies such that the energy 
of the k^^ event is E"^ = Ek± (j{Ek; Psys)- For the systematic errors we choose 
pi = 0.05 and p2 = 0.10. 

To account for this in the parameter estimation procedure, we shift each en- 
ergy up or down and carry out the methods in §3.1[ The difference between 
the parameter estimates of a shifted set and those of the centered set gives 
"systematic" errors of the parameter estimates. 



3.3 Statistical Energy Error 



To model the statistical energy errors of the detector we assume that the 
true energy of the cosmic ray has a 68% chance of being within the interval 
{Eobs — o'stat, Eobs + o-stat)- The observed energy has been "smeared" from the 
true value; Eobs = Etme + ^ where Y is drawn from a normal distribution with 
mean and variance (Tstat- Note that while the true energies can only be found 
on [-Emin, 00), there is a nonzero probability for the (after smearing) observed 
energy to be less than -Emin; Eobs lives on the interval (—00,0x3). This edge 
effect near E^am can be accounted for by assuming that the true distribution 
of energies follows a power-law well below -Emin and then re-normalizing the 
convolution technique used in Howell^. See ^for further discussion. For the 
integrand, three factors are necessary: 

(1) The model to be fitted, fuitj) (see g23D. By letting ^0 = O.l^min we 
are assuming that the power-law extends below the observed -Emin- 

(2) A normal distribution G(t; E^hs, o'statit] P)) with mean E^hs and variance 
<^stat(t'iP) to reflect the statistical energy errors. 

(3) The acceptance of the CR detector as a function of the true energies fl{t). 
Since we are using MC data we choose Q{t) = 1 for simplicity. 



^ With its hybrid detector the Auger reduces the systematic error to between 7% 



and 15%[15j. 
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The convolution is calculated by integrating over all possible true energies it): 

gM{E,^s;0,p)= fu{t]d)G{t-E,^,,astat{t\P))^{t)dt. (3) 

Re-normalizing so that the observed energies define a p.d.f., we numerically 
calculate the p.d.f. to be: 

jM[Eohs,(^,P) - — , , , l4j 

and we must modify the likelihood found in Eq([T]) accordingly: 

CMi9)=Y.ln{fM{E,;9)}. (5) 
1=1 

By finding the parameters 6 which maximize Eq^S]) we can be confident that 
we are accounting for the statistical uncertainty inherent in data collected 
by a realistic detector. To model statistical errors in our toy data set, we 
parameterize astat as in Eq([2]) with pi = 0.15 and p2 = 0. 



4 Evaluating the Fit 

In this section we outline ways to evaluate the fit of a candidate model to the 
data set. The Kolmogorov-Smirnov statistic can be used to extract a best fit 
minimum energy -Emin and, with its corresponding p-value, evaluate the "ab- 
solute goodness of fit" of a candidate model (see §4.ip . The relevant question 
for CR physics is not whether a particular model is a good fit to the data but 
rather whether the fiux exhibits suppression (relative to the single power-law 
form) at the highest energies. To address this question directly we use two 
statistics with well defined p- values: the Tail- Power statistic (see §4.2p . which 
can give information about tail suppression in standard deviations, and a like- 
lihood ratio that allows rejection of the single power-law hypothesis in favor 
of a suppressed candidate model (see §4.3p . 



4-1 Kolmogorov Statistic 

While the minimum value of the likelihood function will indeed give the best 
value of the fit parameters, this fit may nonetheless be poor. The typical [SP] 
method for evaluating goodness of fit is the Kolmogorov-Smirnov test [IS]. The 
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relevant statistic for this test is the KS distance: 

/^Ks(^min) = max |Ffit(E) - Fdata(^)| , (6) 

where, Fat and Fdata are the cumulative distribution functions (c.d.f.) of the 
best fit model and the data respectively. The maximum distance between 
the c.d.f. 's is taken over all energies in the fitted data set, E > -Emin- By 
stepping over -Emin and re-minimizing Eq([T]) at each step to determine the 
best fit parameters, we can calculate Dks as a function of -Emin- The value 
of 00 = -Emin that minimizes Dks can be taken as the best estimate of the 
minimum energy above which the model holds [1]. 

To test how well a particular model fits the data we must simulate many MC 
data sets drawn from the best fit model p.d.f. with the same number of events 
as the original data. The fraction of sets pxs with Z^ks greater than that of 
the data gives the suitable p-value; if pxs ^ 1 then it is unlikely that the data 
are drawn from the model under consideration, and in this way the KS test 
statistic pks can rule out the different candidate modelsflj. 



4.2 Tail Power Statistic 



The Tail-Power (TP) statistic is similar to the KS statistic discussed above, 
however it has, at least, three advantages over pks when testing the power- law 
assumption; 

(1) The TP statistic and it's corresponding p- value ptp are nearly indepen- 
dent of the value of the spectral index 7, 

(2) The asymptotic behavior of the TP statistic is known, and therefore no 
simulations are required to calculate the corresponding p- value pxp, 

(3) If TP > the deviation suggests flux suppression in the tail and if TP < 
the deviation suggests flux enhancement in the tail^ and 

(4) ptp offers an unambiguous p-value in standard deviations. 

This "measure of power-law-ness" has been developed and studied elsewhere 
(see p!2l[T5|7] ) and here we expand its use to the un-binned case. 

The sample TP statistic is defined as [T^ : 

f(£^mm) = Z>i(Emin) - ^ Z>2 (-E^min) , (7) 

where: 

^n(i?mi„) = ^ E (8) 
-^^> E ->E ■ -^min 
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and the sum is carried out over all A^> events with energy greater than a given 
minimum. If the data are drawn from a pure power-law then f{E^[^) will tend 
to zero as — > cx), regardless of the value of 7[5J. 

We may approximate the asymptotic joint distribution of z>i and U2 as a bi- 
variate Gaussian /^^^^(z/i, z/2). The asymptotic mean and variance of ui are 
and and of z/2 are (^zriya and j^^^zij^- The random variables ui and z/2 

are highly correlated; the correlation coefficient is p = independent of^y. 
Thus, for a given and 7, we calculate the p.d.f. of r to be, 

/oo 
U,,{t,2{t'-T))dt. (9) 
-00 



The analytic "location" {t)tp ~ and "shape" {<Jr)TP = \J{t'^)tp — {'t)tp ~ 
jV~i/2(^^ _ l)-2 parameters of this distribution are consistent with simulation 
generated values. We measure the p-value p^p for the TP statistic in units of 
standardized deviation, 

/T-i \ ^(-^min) — (t) TP /-.„N 
PTP(-Emin) = -, ^ • (10) 

{(^t}tp 

A spectrum with flux suppression in the tail (like that in the Fermi-like model) 
will result in a positive significance [7]. 

The application of Eq ffTOl) to the toy CR data set (see §2.21) is plotted in 
Figini The top panel shows the (pure power-law) spectral index as a function 
of -Emin- A spectral index which increases as -Emin increases is indicative of fiux 
suppression. The red, left leaning hatching shows the variation of 7 due to a 
±l(j systematic shift in the energies (see §3.21) while the opposite, blue hatching 
shows the statistical error of the estimator 7, see §3.1[ The bottom panel shows 
the resulting TP statistic significance PTp(-E'min) in standard deviations. Notice 
that while the systematic errors can be significant for the measured spectral 
index, they do not effect the TP statistic. Since we must estimate the spectral 
index to compute ptp, we also propagate the statistical errors on 7 to the tail 
power statistic. 

To test the effectiveness of this statistic, we apply it to a series of simulated 
data sets drawn from both the Fermi and double power-law models. For all 
the models we se0 -Emm = l.OEeV, 7 = 2.75 and either S = 4.75 or Wc = 0.10. 
We vary each characteristic cut-off energy, either E]^ or Ei, in three steps 
lg(ii^cut/-£'min) = 0.5, 1.0, and 1.5. The total number of events in the data set 
is varied in four steps Ig(A^) ~ 2.5, 3.0, 3.5, 4.0. For each of these twelve sets 
of parameter choices we make 10^ Monte-Carlo realizations and plot the mean 
and RMS of pTp(^mm = 1-0) in FigJl 

^ These values are similar to the Auger [15] and HiRes[l] best fit values. 
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Fig. 3. Top The best fit (see ^3.11 Eq([T])) spectral index 7 as a function of Ig E'min for 
the the toy CR data set (see §2.21) fit to the pure power-law model (P). Bottom The 
resulting TP statistic significance PTp(-Emm) in standard deviations as a function 
of the minimum energy i^rniii, see Eq fllOp . Both plots give strong evidence of flux 
suppression of the highest energy MC events. 




Fig. 4. The tail power significance, PTp(-E'min = 1-0) as a function of the (log^g of 
the) number of events in each Monte-Carlo realization. Each plot style represents 
a different choice of lg(-E'cut/-E'mm) = 0.5, 1.0, or 1.5. Le/t, the double power-law, 
Ecut = ^b- Right, the Fermi power-law, i^cut = Ei. 



Based on FigJH we can see that the best way to evaluate a data set with a 
potential for tail suppression is to collect as much data with £"111111 as close to the 
expected cut-off as possible. The experimenter may use FigHJ or one like it, to 
help tune observation parameters, i.e. collecting time on a gamma ray source 
or size of a CR detector, in advance of the observation and in anticipation 
of flux suppression of a certain type. Note, however, that one should choose 



10 



an Eynin prior to analyzing a data set to avoid a penalty for scanning in this 
parameter. 



4.3 Model Discrimination 



Here we introduce a likelihood ratio test designed to discriminate candidate 
suppressed models (DP and FP) from the pure power-law. We define two 
log-likelihood ratios; for each model M: 



N 



7^M = E {^MiEi) - £p(E,)} = Cm- Cp, 



'111 



i=l 



where iuiEi) = \iafM{Ei;9) with M either DP (double power-law) or FP 
(Fermi-like), and ip{Ei) = lnfp{Ei;6) for the pure power-law likelihood per 
event (see Table [1] and Eq([T])). Note that each suppressed model is fit inde- 
pendently of the pure power-law best fit. The asymptotic variance of 71 can 
be estimated by the sample value: 



N 



N 



Y,\[iMiE,)-ip{E,)] 



N . 



(12) 



The hypothesis of the pure power-law is nested within the hypothesis of a 
suppressed power-law. As a consequence, |7^|/cr7^ ^ 0/0 as — 00 and the 
distribution of TZ/a-n is not Gaussianjl]. The correct p- value is calculated as 
the integral of a function |14|H] : 

Pniz') = ^ r t~'/'e~'/'dt, (13) 

where = 7^^/(2A^a|). 

We interpret this p-value in the following way: if p-ji is "small" then the best 
fit model M may be preferred over the best fit pure power-law. By small we 
mean that, a priori and rather arbitrarily, we may choose to reject the single 
power-law in favor of the model if pn < 10^^. This quantity tells us only 
whether a given suppressed model is better than the pure power-law. It says 
nothing about how well any of the fits actually represent the data. 

For each of the twelve sets of parameter choices used in FigSJ we plot the mean 
and RMS of pn in FigS As before, we see that the best way to reject the 
power-law in favor of the suppressed model is to collect as much data with -Emin 
as close to the expected cut-off as possible. Note that for \g{Ecut/ Emm) = 1-5 
the distribution of likelihood ratios is strictly positive and highly peaked near 
zero; the mean and RMS are not good reflections of this distribution. 
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Fig. 5. The log of the hkehhood ratio significance, p-ji as a function of the (log^^o of 
the) number of events in each Monte-Carlo realization. Each plot style represents 
a different choice of lg(-E'cut/-E'min) = 0.5, 1.0, or 1.5. Left, the double power-law, 
£^cut = E\j- Right, the Fermi power-law, Ecut 
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5 Summary and Conclusion 



In this paper we describe a set of statistical tools designed to extract the 
most accurate and precise information about the flux of the highest energy 
cosmic rays. We show how to use the un-binned likelihood method described 
in §3.11 to fit a data set to the three model distributions described in §2.31 
Techniques for incorporating the systematic and statistical errors associated 
with a real CR detector into the likelihood method are described in §3.21 and 
§3.31 respectively. In §1] we describe p- values useful for extracting information 
about flux suppression. We show in §4.21 and §4.31 how an experimenter might 
use an a priori estimate of the cut-off energy to maximize an observational 
setup for detecting flux suppression. 

The collection of these statistical tools are the primary result of this paper. To 
answer the questions posed in the introduction for a given data set we suggest 
the following steps: 

(1) Estimate the best flt parameters 6 of the model; 

(a) The estimates 7, E'b or Ei and 6 or lUc are determined via the like- 
lihood Eq©, 

(b) The estimate of the minimum energy i?min is that which minimizes 
the Kolmogorov distance Dks (see §4.11) . 

(2) Shift the energies up and down according to the systematic uncertainty 
described in §3.21 and repeat step ([1]) . The resulting shift in parameter 
estimates gives the systematic uncertainty of those estimates. 

(3) Obtain the model parameter estimates using the methods in §3.31 to in- 
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corporate the statistical error of each event energy. 
(4) Test the model hypothesis; 

(a) The absolute goodness of fit for any of the models can be evaluated 
using pks in §0, 

(b) The Tail-Power statistic ptp can be used to reject the single power- 
law hypothesis (nearly independently of the spectral index estimate, 
see 

(c) The single power-law may be rejected in favor of a specific alterna- 
tive model using pn, here we study the double and Fermi power- law 
distributions (see §4.31) . 

The best estimates for the characteristic cut-off energy and shape parameters, 
determined via steps ([1]), ([2]) and ([3]), are £'b or Ei and 6 or lUc respectively. 
The presence of flux suppression at the highest energies can be evaluated using 
step dH). 

By applying these methods to the toy Monte-Carlo set of CRPropa events 
we illustrate in ^ how the procedure may be implemented on an actual 
CR detector, i.e. a detector with systematic and statistical event energies. 
Suppression in the tail is clear in Fig JC.ll and Fig JC.2j the tail power statistic 
is 4.6(7 and the p-value for the double (Fermi) power-law is Igpop = —2.7 
(IgPFP = -1.9). 

The methods are sufficient and robust. Indeed, many of them have been ap- 
plied by the Auger collaboration which reports suppression with 6a confidence |T3]. 
These tools serve as a basis for further investigation of the CR spectrum such 
as evidence for more detailed spectral information. They can be applied to 
any data set, astrophysical or otherwise, to provide information both about 
data already collected and help to optimize future observations for detecting 
tail suppression. 
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A Binned vs. Un-Binned 



The statistical superiority of an un-binned maximum likelihood estimate of the 
pure power-law spectral index to the logarithmically binned least-x^ method 
often used has been established in [S] and expanded upon more recently in 
[9][Tl|il[71l8] ■ In this section we compare the binned to the un-binned fitting 
method for the two suppressed models, i.e. the double and Fermi power-laws 
(see g23D. 



To calculate the binned estimators we minimize a x^(^) function that relates 
the logarithmically binned (width w) histogram of the data to that expected 
by a model. The function igfj, 

where Nb is the number of bins, Yf^^^ is the number of events in the z*^ bin 
bi and Uj is determined by Gaussian errors when y.'^^^''^ > 10 and Poissonian 
errors when F^'i^*^ < 10. We minimize with respect to the parameters 6 (with 
9q = -Emin fixed) using the number of events in a bin expected by the model 
M, 

Y^\e)=N fM{t;e)dt. 

JlQbi-w/2 

To study the asymptotic bias and error produced by the two estimation tech- 
niques we draw 10^ sets of 5 x 10^ events from a pure power-law and separately 
from a double distribution. For each Monte-Carlo set we estimate the best fit 
model parameters 6 using both the likelihood Eq([T]) and the Eq (]A.ip meth- 
ods. The un-binned estimator of the pure power-law spectral index (see §3.11) 
has been shown [5ll9] to have an error estimate within ~ 1% of the Cramer-Rao 
lower bound for a sample with as few as ~ 100 events. 

In Fig lA.ll and Fig lA.2l we plot the results of the simulations. We can conclude 
that the un-binned fitting method is most important when fitting a power-law 
in the tail of a distribution; the binned estimator performs nearly as well as 
the un-binned for the double power-law parameters 7 and E\y. The (binned) 
methods used to report parameters like the "ankle" and the "knee" in [15] 
and |Tj are sufficient but limited by the bin width. 



^ For the case of the single power-law Ig /p = Ig C — 7 Ig where C is the normal- 
ization. Thus the binned fitting method reduces to fitting the logio of the (error 
weighted) bin heights to a straight line with slope 7. This technique is often used to 
mitigate the effects of the heaviness of the power-law tail but un-binned methods 
are more accurate and precise. 
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Fig. A.l. For each of 10^ sets of 5 x 10'^ events drawn from a pure power-law 
with index £"111111 = 1-0 and 7 = 2.75 we estimate the spectral index using the 
binned Eq ijA.lh and un-binned Eq([T]) methods. The bias and error of the un-binned 
estimator is 0.0002 and 0.0247 and that of the binned is -0.024 and 0.0272. 




T Eb 5 



Fig. A. 2. For each of 10^ sets of 5 x 10^ events drawn from a double power-law with 
parameters {'y,E\i,6} = {2.75,10.0,4.5} we estimate the spectral index using the 
binned Eq ijA.lh and un-binned Eq([T]) methods. The bias and error of the un-binned 
estimators are {—0.002, 0.13, 0.16} and {0.03, 1.4, 0.60} and those of the binned are 
{-0.005,-0.71,-0.43} and {0.03,1.7,0.60}. 
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B Statistical Error: Monte-Carlo Example 

To illustrate the effect tlie statistical energy smearing has on a pure power-law 
we generate 9000 MC events from a power-law distribution with -Emin = 1-0 
and 7 = 2.75. A histogram of these events is represented by the black filled 
circles plotted in Fig lB.li By minimizing Eq([T]), we calculate the estimated 
spectral index for this data to be 7 = 2.742 ± 0.019 (with -Emm = 1-0, see 
§4.ip . A power- law with these parameters is plotted as the dashed line in 

FigEH 

To each MC event Ei we then add a random number Yi drawn from a normal 
distribution with mean zero and variance 0.2Ei. The new events are histogram- 
ed with blue open circles in Fig lB.li We fit these events by maximizing a 
likelihood with 

poo _, 

/ fM{t]e)G{t-E,^,,(Jstat{t\p))dt. (B.l) 

(compare with Eq([3])) as the p.d.f. and we find that 7 = 2.749 ± 0.020. The 
smearing does not effect the estimated spectral index, though it does increase 
the error of the estimate. The dashed curve in Fig JB.ll shows Eq (lB.ip evalu- 
ated at the best fit values. Notice that the histogram of the smeared energies 
deviates from the un-smeared case near Igi? ~ 0. In §3.31 we account for this 
edge effect at the low energy end by assuming that the true energies follow the 
power-law well below the observed minimum energy; in constructing the like- 
lihood we choose O.lE^mm for the lower rage of integration (compare Eq flB.ll) 
with Eq([n])) and we re-normalize to ensure a true p.d.f. (see EqQ). 



C Results of CRPropa Toy Set 

By applying the statistical tools presented in this paper (summarized by steps 
(Pl-dlQ in g5]) to the toy set of 5 x 10^ CRPropa events (see g22D we il- 
lustrate how the tools might be implemented on an actual CR detector. By 
construction, this toy set has parameter estimates and, more importantly, er- 
rors estimates and hypothesis test p-values that are numerically comparable 
with those reported by Auger [15] and HiResjl]. 

In preparation for this paper we generated 14 CRPropa simulations of ~ 2 x 10^ 
events with different injection spectral indexes, 71N = (2.0, 2.1, . . . , 2.6), and 
with different values of maximum generation energy, E'max/EeV = (400, 2000). 
The (after propagation) estimated characteristic break point energy, i.e. Ei or 



Note that since we are not interested in the absolute goodness of fit for any of 
these toy models to this toy data set, we do not perform step (flaj) of ^ 
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Fig. B.l. An example of a pure power-law before and after smearing. A histogram 
of 9000 events drawn from a single power-law with -Bmin = l-OEeV and 7 = 2.75 is 
plotted in black filled circles. The best fit (using Eq([T|)) power-law for these events 
is plotted in solid black. The blue open circles are a histogram of these events after 
being smeared by a Gaussian with variance 0.2E (see ^B]) . The blue dashed curve 
shows the best fit using Eq ljB.ip . To account for the edge effect near Igii^ ~ we 
use the methods in §3.31 namely Eq([3]). 



El,, is found to be independent of the spectral index at the site of generation, 
7iN. The estimated spectral index 7out is found to be linearly related to the 
input spectral index 71N with linear slope ~ 1. The high energy estimated 
shape parameters, 6 and Wc, are more sensitive to the maximum generation 
energy (at the sources) than they are to 7in. 

In Figs. IC.ll and IC.2I we plot the toy data set and the best fit models in two 
(non-binned) ways not commonly seen in the CR literature. The first is a 
rank- frequency plot. For each event (black filled circle) we plot IgE along the 
horizontal axis and the log of the number of events with energy greater than 
E along the vertical. For each of the models (see §2.31] , the vertical axis is 
\g{Ntot{^ — F{E))) where F{E) is the model cumulative distribution function. 
From the rank-frequency plot we derive an instructive visualization tool in 
Fig JC.2j we plot the difference between the number of events above a given 
energy for the toy set NS}^^ and that expected by the best fit models A^>^^. 

The best fit pure power-law parameters for the toy set described in §2.21 are 
^min = 6.31±0±[j;| and 7 = 2.83±[]:[]|±;o:?o where the first error is statistical 
and the second systematic. The tail power significance pxp is 4.6cr. The best 
fit double power-law parameters for the toy set are -Emin = 6.31 ± ± 0.82, 
7 = 2.71±0.03±;Ho, = 45. 7 ±1? ±9.9 and 6 = 4.30±0.26±;[JiJ. The cor- 
relation coefficients are p^E^ = 0.18, p^s = —0.15 and pe^^s = 0.32, see Fig |C.3[ 
The likelihood ratio significance is IgpT^ = —2.7. The best fit Fermi power-law 
parameters for the toy set are E^^^ = 6.31 ± ± 0.82, 7 = 2.69 ± 0.03±;{]:|]^, 
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El = 78.6 ±\li and = 0.139 ±°:°24 ±0-005_ ^^^g correlation coefficients 
are p^Ei = 0.61, p^^^,^ and Peiwc = —0.07, see Fig JC.4[ The likelihood ratio 
significance is IgpT^ = —1.9. 




lg(E . /EeV) 

mm ' 

Fig. C.l. A rank-frequency plot as simulated by 5 x 10^ events from the CRPropa 
set with parameters 71N = 2.6 and i?max = 2000 EeV. For each event (black filled 
circle) we plot Ig E along the horizontal axis and log-number of events with energy 
greater than E along the vertical. The models are described in ^2.31 
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Fig. C.2. Using the rank-frequency plot (see Fig lC.ip we plot the difference between 
the number of events above a given energy for the toy set A^"*^*^ and that expected 
by the best fit models N^^. Note that at lg£'min/EeV ~ 1.7, there are at least forty 
fewer events observed than expected by the pure power-law fit, i.e. flux suppression. 
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Fig. C.3. The change in log-likehhood — 2A£dp (see ^3.ip as a function of the pa- 
rameters 7, -Eb and 5 of the double power-law. The data set is the toy set described 
in ^2.2[ The best estimate for each parameter is plotted as a blue box, the asym- 
metric one degree of freedom error estimates (— 2A£dp = 1) are plotted as solid 
blue lines and the black contour defines the two degree of freedom error estimate 
(-2A£dp > 2.30). 




Fig. C.4. The change in log-likelihood — 2A£fp (see §3.ip as a function of the 
parameters 7, Ei and Wc of the Fermi power-law. The data set is the toy set 
described in ^2.2f The best estimate for each parameter is plotted as a blue box, 
the asymmetric one degree of freedom error estimates (— 2A£fp = 1) are plotted 
as solid blue lines and the black contour defines the two degree of freedom error 
estimate (-2A£fp > 2.30). 
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